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Abstract. We present a description for setting initial particle displacements and field values 
for simulations of arbitrary metric theories of gravity, for perfect and imperfect fluids with arbitrary 
characteristics. We extend the Zel’dovich Approximation to nontrivial theories of gravity, and show 
how scale dependence implies curved particle paths, even in the entirely linear regime of perturba¬ 
tions. For a viable choice of Effective Field Theory of Modified Gravity, initial conditions set at high 
redshifts are affected at the level of up to 5% at Mpc scales, which exemplifies the importance of going 
beyond A-Cold Dark Matter initial conditions for modifications of gravity outside of the quasi-static 
approximation. In addition, we show initial conditions for a simulation where a scalar modification 
of gravity is modelled in a Lagrangian particle-like description. Our description paves the way for 
simulations and mock galaxy catalogs under theories of gravity beyond the standard model, crucial 
for progress towards precision tests of gravity and cosmology. 


Contents 


1 Introduction 1 

2 Density displacement duality 3 

3 Linear displacement and velocity fields in metric theories of Gravity 5 

3.1 Geometry 5 

3.2 Hydrodynamics 5 

3.3 Constant or variable mass per vertex during the fluid phase 6 

3.4 Bare densities and velocities 7 

3.5 Adiabatic and isocurvature initial conditions for multiple fluids and fields 8 

4 Non-relativistic initial conditions 9 

4.1 Non-relativistic limit 9 

4.2 Gauges and the Newtonian limit 10 

5 Scale dependent growth in the quasi-static approximation 11 

6 Effective Field Theory parametrisation, and modified gravity as a particle dis¬ 
placement field 13 

6.1 Transfer functions 13 

6.2 Initial conditions for discrete simulations 18 

7 Conclusion 19 

A Gauge transformation and Q-fiuid variables in the longitudinal gauge 20 

B Discrete sampling of phase space 21 

B.l Finite length 21 

B.2 Finite resolution 21 

B.3 Generation in Fourier space 22 

B.l Generation in real space 22 

B.5 Pre-initial conditions, or particles 23 

B.6 Summary 23 


1 Introduction 

The expansion history of the universe seems to be a smooth function of time, which can hence only 
provide a handle on a limited number of parameters in the cosmological model. In order to further 
quantify Inflation, Dark Matter and Dark Energy, and distinguish a variety of candidate models, hope 
is placed on probing the nonlinear growth of perturbations in the universe, considered the ongoing 
preparations for the Euclid satellite [1]. 

The current state of the art for nonlinear structure formation in the universe is the numerical 
simulation of gravitational clustering of a large fixed number of particles in a finite sized box: N- 
body simulations. A plethora of codes solves Newtonian dynamics endowed with a friction term that 
accounts for the cosmic expansion {e.g. [2-4]). Progress is made in relativistic simulations as well [5]. 
The next layer of complexity is the addition of perturbations in the exotic form of matter mentioned 
above, Dark Energy, or often under the alternative moniker of Modified Gravity [6-10]. 

In order to push simulations to the next level, accounting for dynamical modifications of gravity 
or Dark Energy at arbitrary redshifts, an obvious first step is to discuss the starting point of such 
simulations, which is the goal of this paper. 
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Numerical simulations are inevitably discrete in space and time, and it is not clear what the con¬ 
sequences of this discreteness are for the results, or how closely the results reflect nonlinear structure 
in a continuous universe [11-18]. Such discreteness effects are not the subject of this paper. However, 
in Appendix B, we summarise the known methods for realising discrete samples of continuous fields, 
and we list the choices that the simulator needs to make. 

In most of the literature on simulations, priority is given to the details and technicalities of 
solving for the nonlinear dynamics. Setting the initial conditions is always done using Zel’dovich’ 
approximation (explained in the main text), sometimes up to second order [19]. Most importantly, to 
our knowledge, all simulations set their initial conditions at a time when no perturbations in exotic 
matter are present. However, if one wants a more complete description, the full parameter range of the 
exotic models needs to be tested, including the range in which Dark Energy has significant dynamical 
perturbations (see [20] for a review); where the quasi-static approximation for its perturbations breaks 
down [21-24]. 

The aim of this paper is to describe all that is necessary to set the correct initial conditions for 
arbitrary cosmological simulations under arbitrary metric theories of gravity, varying from the inclu¬ 
sion of only dust to the inclusion of species that in their linear description are fully imperfect fluids. 
For the sake of clarity, we include discussions of previously known topics where necessary, providing 
the relevant references. We do not address the subject of setting initial conditions for matter species 
that do not simply form a sheet in 6-dimensional phase space, such as for example standard-model 
neutrinos and photons above a certain resolution, which are described accurately by linear perturba¬ 
tions although not by fluid dynamics. We release a code, FalconIC which can be used to generate 
discrete initial conditions of any cosmological fluid. The code links against any version of both Boltz¬ 
mann codes CAMB [25] and CLASS [26], including EFTCamb [27, 28], such that no separate running 
of those codes is necessary, and initial conditions for any fluid can be generated at arbitrary scales, 
of arbitrary size (fully parallelised using MPI and OpenMP), for arbitrary cosmological parameters. 

Simulations discretise the cosmological fluids in a regime where the fluid description is still valid, 
on the onset of shell crossing and the onset of the need for a particle description, when the sheet in 
phase space starts folding or even tearing. Therefore, at the first time steps, N-body simulations should 
reproduce the linear theory described by fluid dynamics. However, the fact that the linear discussion 
in this paper addresses (imperfect) fluids, does not at all imply that these discrete realisations can 
only be used for fluids. In other words, the Effective Field Theory of Modified Gravity describes some 
matter which in the linear regime may look like an imperfect fluid. Particles whose free-streaming 
length is smaller than the simulation resolution, but whose dynamics is described as in imperfect fluid, 
such as standard model neutrinos, can be realised just as well following our description. 

In summary, the following points are the ingredients for discrete realisations of arbitrary imper¬ 
fect fluids, 

• Any linearly perturbed quantity (‘charge’) can be translated into a displacement field of equal 
charge vertices, by a coordinate transformation whose Jacobian equals the original charge per¬ 
turbations. 

• The scalar quantity that defines the positions of particles with velocity for an arbitrary fluid, 
is = ~ f drWiu'^, where i runs over spatial coordinates. 

• For newtonian simulations, positions are given by n = /\fg^, with the determinant 

of the spatial part of the metric. 

• In the presence of pressure or heat flux, the particles have varying masses m and internal energies 
r, defined by 

m (I + Ap - A„) = ^ (1 + A^“^ - A]]“^) , (I.l) 

n n ^ ' 

T =- (1 + Ap - A„) = - (1 + - A^”'’) , (1.2) 

n n 

^http://falconb.org 
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Figure 1. Approximating Zel’dovich: a discretised displacement field based on an image of Y. Zel’dovich plus 
gaussian noise, multiplying the displacement with increasing factors from left to right, starting at zero. 


with Aq = 6a/a for any quantity a with background value a. 

• For adiabatic perturbations, all quantities are set by the same random seed multiplied by their 
respective transfer functions. Isocurvature perturbations are introduced by combining multiple 
random seeds. 

• The linear displacement field for a dust fluid in a universe endowed with a scalar modification 
of gravity is scale dependent, and hence the trajectories of particles are not straight lines. 

• The scalar modification itself can be described in both Lagrangian and Eulerian representations, 
whichever of the two may prove more convenient in nonlinear simulations. 

The bending of Dark Matter trajectories at the linear level, implies that simulation codes such 
as Ramses [2] or Enzo [4] which only take particle velocities as input, need to be adapted to include 
particle positions as well, since particle positions are no longer simply a time-dependent function times 
their velocities. 

In section 2 we start by explaining how an arbitrary density field is related to a displacement 
field, under the assumption of absence of vorticity^. In section 3 we then show how to define the 
displacement field under an arbitrary metric for an arbitrary (im)perfect fluid, with the correct ve¬ 
locities emerging. In section 4, we briefly comment on the Newtonian approximation, necessary for 
Newtonian simulations. In section 5 we show how particle trajectories of linear perturbations do not 
follow straight lines, already in the quasi-static approximation. Finally, in section 6 we apply our 
prescription to a parameterisation of the Effective Field Theory of Modified Gravity, in synchronous 
comoving coordinates. 

Throughout this paper, we will refer to the Conformal Newtonian gauge as the Longitudinal 
gauge, since we prefer to reserve the word Newtonian for Newtonian N-body simulations.We use units 
with the speed of light c = 1, and we use the Einstein summation convention. All equations with 
perturbative quantities are expanded up to linear order. 


2 Density—displacement duality 


A displacement field can be regarded as coordinate transformation from some (virtual) coordinate 
system with a constant charge per coordinate volume to a (physical) coordinate system with a per¬ 
turbed density (of charges, mass, etc.). When the perturbations are small, and the unperturbed 
quantity is denoted by an overbar, the coordinate transformation to obtain the displacement field 
that is associated with some density field p{x) = p{l + D{x)) can be defined by. 


p{l + D{x))d^ ^x=pd^ ^x', 


{l + D{x)) 


dx^ 
dx^ ’ 


^ One can add vorticity by means of several scalar fields [29] . 


( 2 . 1 ) 

( 2 . 2 ) 
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such that an observer at fixed x' is co-moving with the charges p. For small D{x^) and in absence of 
vorticity, the coordinate transformation is solved by 


where = Yli' ^1' 


X* =x*' - :^D{x') 

=x^' -^D{x)+0{D^), 


1 


D{x) = 


^x 



ik-{jc—x^ 


(2.3) 


(2.4) 


This relation is independent of the theory that gives rise to the density field^, as exemplified in Fig. 1. 
However, since we set vorticity to zero, the Poisson equation that determines the Newtonian potential 
is equal to the equation that gives the deformation tensor dx'' jdx^ . One can identify the scalar D(^x) 
with the Newtonian potential, and obtain the Zel’dovich Approximation [30, 31]. 

In the x' coordinate system, the equal-weight particles are not displaced, such that this is a 
Lagrangian coordinate system^: the coordinates are at rest with the charges, which do not necessarily 
have to be masses. 

If the density field has a continuous time dependence, p{x) p{t,x), the velocities of the 
virtual equal-weight particles are simply given by the time derivative of the field, 


#‘=)(T,y) = -^dtp{T,x). (2.5) 

When one is dealing with a curved manifold, care needs to be taken with the meaning of the 
displacement and density fields. For the displacements to correspond to the coordinate positions of 
particles on an arbitrary manifold, the density field used for the coordinate transformation must be 
the ‘bare’ density of particles in coordinate space, not taking into account the curvature of space. If 
one were to take the proper density perturbations as the generator for the coordinate transformation, 
one would obtain the displacement of particles in a (purely spatially) transformed coordinate gauge, 
in which the trace of the perturbations of the spatial metric vanishes. 

Instead of employing the purely spatial density-displacement duality above, one may perform a 
similar coordinate transformation using the proper energy density on the spacetime, 

p{x)V^d^x = p{l + D{x))y/^d^ X = p\/—g'd^x', (2.6) 


in which case the x' coordinate system acquires the meaning of a synchronous comoving gauge, 
associated with the species of p, as explained in [32-34] . Such a transformation makes explicit where 
General Relativity comes into play compared to Newtonian dynamics, but is not useful for the purpose 
of this work. 

When a problem is discretised (see Appendix B), the Eulerian description amounts to following 
a density field on a regular grid in x at positions Xijk, while the Lagrangian description amounts to 
following a displacement field for points on a regular grid in x', which translates to a curved mesh 
in x-space. Here ‘regular’ does not necessarily mean Cartesian, because for example glass initial 
conditions [35] or alternatives [36] take an irregular partitioning of x'-space. 

In general, at any time dM{T,x) = pd^x' = p d^x. This expression continues to hold in 


situations where x'(x) is not single valued and S{x) > 0{1), i.e. the phase-space sheet has folded [37]. 
It must be noted that its evaluation becomes nontrivial then. 


similar conclusion is drawn in Ref. [13], based on the continuity equation such that velocities are in agreement 
with a physical theory. 

^In cosmology, the labels x' are often referred to as the initial positions of particles at the past singularity. This is 
not a consequence of Lagrangian perturbation theory, but this is a consequence of the fact that in cosmology usually 
only growing modes are considered, such that D{x) —> 0 at past infinity. Fundamentally, it is not necessary for x’ and 
X to coincide at any time, especially since in the context here there is no reference to time at all. 
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3 Linear displacement and velocity fields in metric theories of Gravity 

3.1 Geometry 

The next section applies to any metric theory of gravity, with a covariant derivative defined by, 

vx=V + ^>^ (3.1) 

{dfigxa + dxgatJ. - dxg^a ), (3.2) 

with the Christoffel connection We do not refer to the equations that source the metric with a 
matter configuration. That is, we define all modifications of gravity as those that change the Einstein 
equations and / or those that can be written as additional matter content inside the energy-momentum 
tensor. 

3.2 Hydrodynamics 

At the linear level and hence at the Canchy surface of N-body simulations, the cosmic density fields 
under consideration can be described as fluids. Let us hence first layout the definitions for relativistic 
hydrodynamics. All contents of this section can be traced back to Refs. [38-41]. 

This paper focusses on dynamics entirely attributable to scalar quantities. All is perturbed 
about a background Friedmann-Lemaitre-Tolman-Bondi solution, denoted by overbars, such that 
background quantities have only time dependence while perturbed quantities (indicated by A) have 
time and space dependence. Hereafter we drop the time and space dependence in most functions. For 
a fluid with four-velocity = dx^/x/—ds^^ U^Ufj, = —1, define the transverse projector®, 

-Lfj.iy= gfiv + (3.3) 

which projects into the plane orthogonal to C/^, the spatial slices for an observer comoving with the 
fluid. The unperturbed (7* = 0, and 9^(517* = 0. 

The energy momentum tensor then carries the following information, 

• energy density p = p + 5p = p{l + Ap) = 

• pressure P = P + SP = P{1 + Ap) = ^ 

• energy flow (or heat transfer) U^T^x, 

• anisotropic shear perturbation 
and is decomposed as 


=pU^^U, + P + U,q^ + E'^,. (3.4) 

Owing to the scalar nature of the system one can further simplify with 

q>^=-a{T)ip + P)A^’^^q, (3.5) 

= - |a(T)2 (p + P) {±^x±.c -I (3.6) 

where we chose pre-factors for convenience, such that a corresponds to the definition in Ma & 
Bertschinger [40] when q = 0. 

In an arbitrary gauge, the scalar part of the metric can be written as [42] , 

, = —(1 -I- 2A)dr^ — 2BidTdx'^ + [(1 + 2HL)r]ij + 2/i^] dx^'dx^ (3.7) 

ayr) 

®The symbol ± can be pronounced ‘perp’, for ‘perpendicular’. 
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where r]^^, is the Minkowsky metric, and 


B, = 


Kj = 


ik-x 


(27r)" k ^ 


^ _ 1 

V2 3 


V^3 


Ht, 


(3.8) 

(3.9) 


where all scalar potentials A, B, and Ht are small and spacetime dependent®. 

Any number density in the frame of an observer comoving with the fluid at velocity is a 
scalar n = n + 6n = h(l + A„), and the number transport is = —nU^. If the number is conserved, 
we have at the linear level and independent of whether the fluid’s energy is conserved. 


V^n'^=0, [number conservation] (3.10) 

h oca“®, (3.11) 

A„ = — (0 + SHl), [linear order] (3.12) 


where an overdot denotes a derivative with respect to conformal time r. The number can be associated 
with microscopic particles, but just as well with macroscopic simulation vertices (also often referred 
to as ‘particles’). 

The energy conservation equation V = 0 corresponds to the continuity and Euler equations, 
at linear level in Fourier space' 


P = -m (p + P), 


(3.13) 

.p = (! + «;) (q-e-mL 

)+3H{w- cl) Ap, 

(3.14) 

q =k^A + k^Ap — H(1 — 

3w)(9 + kB — q) - (0 + kB — q) — k^a, 

w -\-l 

(3.15) 


where H = d(T)/a(r) and w { t ) = Note that at the linear level, the continuity equation is only 
sensitive to the trace of the spatial perturbations of the metric, 3Hl- The scalar potentials A and B 
only affect the Euler equation for the velocities, while the transverse traceless potential Ht does not 
enter the energy conservation equations at the linear level. 


3.3 Constant or variable mass per vertex during the fluid phase 

The phase-space of a fluid is discretised for the purpose of a numerical simulation, upon the start of 
which a particle picture may be followed, which goes beyond the fluid description when perturbations 
become nonlinear and trajectories start crossing. This means that a mesh is laid out in real space 
with a velocity associated to each vertex of the mesh, the ‘particles’. For a summary of the existing 
methods for discretisation, see Appendix B. 

Although this it is not a fundamental condition, all numerical cosmological simulations known 
to the authors are based on a fixed number of vertices in phase-space (particles) which is preserved 
throughout the simulation. This is not to be confused with the Poisson solver, which during the 
simulation at any time can choose an adaptive mesh to approximate the gravitational potential(s) at 
the desired resolution, provided the positions of the fixed number of particles. 

If the vertices of the fluid are comoving with the velocity field of the fluid, then their proper 
number density follows Eq. (3.12). Comparing Eq. (3.12) to Eq. (3.14), it is clear that the special 
case of a fluid with constant equation of state ri; = 0, constant sound speed w = which for this case 
is equal to = SP/6p, and zero energy flow g = 0, can be modelled discretely by only accounting for 
the positions and velocities of vertices, associating an energy density to each vertex with, 

Ap = (1 + w)An. [w = q = 0,w = cl = 0 ] (3.16) 

^Compare with the convention in [40], in the longitudinal gauge we have (A = 'i/;, = —4>) and in the synchronous 

gauge = h, Ht = 6ri h). 

^See Ref. [43] for the special case of the perfect fluid, i.e. q = a = 0. 
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Note that the fluid need not be perfect, as tr need not be zero for this condition, which is why the 
statement above may apply to other fluids than dust. The quantity Ap is to be evaluated on the 
discrete vertices on which initial conditions are generated. 

If any of the three conditions is broken during the time span of the simulation, one needs to model 
the dynamical energy density of the vertices. In other words, the ‘particles’ do not have a ‘constant 
mass’, and energy density and pressure at each vertex are given by equating p(t,x) = m{T,x)n{T,x) 
and p{t,x) = n{T,x)T{T,x), where m can be thought of as mass and T as temperature, although 
strictly speaking they can have other meanings, 

m=|(l + Ap-A„), (3.17) 

P 0 f 5P \ 

T =— (1 + Ap - A„) = 3 ( w + — Ap - wA„ ) . (3.18) 

n n \ op J 

It should be clear now that m and T are evaluated at the discrete set of vertices (‘particles’). One 

is free to set for example ^ = 1 at some given time. Moreover, depending on the extent to which a 
gauge is fixed, one may have further freedom to set A^ = 0 at a convenient time of choice, which 
amounts to different choices of cutting the density field into particles. 

3.4 Bare densities and velocities 

A simulation acts in a coordinate space. The bare density of vertices is at any time step given by 
the number of vertices inside a given volume in coordinate space [44]. This is related to the proper 
density by, 


^ n^/gi^d^X, (3.19) 

where = det pij = a(T)®(l + GiJp) + 0{e^), the determinant of the spatial part of the metric, the 
three-metric, such that® 


A„ = A^-- - 3Hl. (3.20) 

Thus, the displacement fields are expressed in terms of coordinates (and not proper distances), such 
that they must be generated in coordinate space based on the bare densities. Density fields, on the 
other hand, continue to express proper densities and need no notion of bare density. 

A displacement field based on number density n, must be consistent with the theory at any time 
step. This implies that the motion along a displacement field at different time steps must reproduce 
the correct velocities. Indeed, we find, 

6 »bare = = -A„ - 3Hl = 0, (3.21) 


as per Eq. (2.5) and Eq. (3.12). 

Note that the definitions of mass and temperature are unaffected by the notion of bare density, 

m = ^ (1 + Ap - A„) = ^ (1 + Aba- _ ^bare^ ^ (3 33 ) 

n n ^ 

T =- (1 + Ap - A„) = - (1 + Aba- _ ^bare) ^ (3 33 ) 

n n 

In summary, particle positions in a relativistic simulation are obtained from applying the density- 
distance duality to bare number densities, while particles masses and internal energies are set through 
the equations above. Again, these quantities are to be evaluated on the discrete positions of the 
simulation particles. 

®As noted in [44], for simulations in the longitudinal gauge, corresponds to / yj Sffomov) 

the comoving gauge, but this coincidence does not occur for other gauges in which the transverse traceless perturbation 
of the metric does not vanish. 
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3.5 Adiabatic and isocurvature initial conditions for multiple fluids and fields 

Now that we have shown in brief how to relate a solution to the stress-energy tensor to a displacement 
field for a quantity described in a Lagrangian picture, in arbitrary coordinates, we can address the 
problem of discretising Lagrangian quantities and Eulerian quantities simultaneously, even though 
the perturbations are strictly at different coordinates. The coordinates {t, x\ are the coordinates 
of the simulation. An Eulerian quantity’s phase space is sampled and simulated on a regular grid 
in T-space, i.e. on Xijk, with indices labeling discrete positions in dimensions a;^, 

respectively. Tracers of the phase space of a Lagrangian quantity are on a regular grid in that 
quantity’s homogeneous-coordinate-density system (“equal-charge tracers”), yijk, with indices {i,j, k} 
labeling discrete positions in dimensions The coordinates are related hy y = x + K, such 

that a regular grid in one space becomes a curved mesh in the other; the grid points in both spaces 
do not coincide. 

How does one generate realisations of initial conditions in real space, when multiple grids are 
present, as many as there are Lagrangian quantities plus one for all Eulerian quantities? The short 
answer is: the same random numbers can be used on all grids, as long as displacements are sufficiently 
small. In a more General Relativistic language: as seeds of linear perturbations, the random numbers 
are gauge independent in Appendix B). This follows straightforwardly from the argument in 
section 2, which can be summarized by Xijk = Hijk + = Hijk + ^{xijk) + where e 

denotes all quantities that are assumed small, being displacements, potentials, densities and so on. 

By virtue of the closed set of linear differential equations that describes the system in the linear 
regime, 

M(T,fc)%(r) = 0, (3.24) 


where M ^r, fcj denotes a matrix with differential operators in t and nonlinear functions in fc, and fj:('r) 
denotes the set of dynamic degrees of freedom, all random input in a realisation of initial conditions 
is encoded in the initial conditions of the differential equations. That is, the differential equations are 
deterministic, and the solution changes linearly with the initial conditions. 




^ o. (-.«)/3 ^ 


(3.25) 


where D ^r, solves the system of equations (3.24). Generating one realisation at any time in cosmic 

history, amounts to choosing random numbers ft"’ that in the real universe actually are drawn at the 
hot big bang, for example by the infiaton as quantum correlations decohere and become classical. 
The quantity ft"’ can be identified with ^ in Appendix B. 

A realisation is entirely described by its Fourier transform, regardless of whether the random 
numbers are drawn in Fourier space or in real space (as in [45] , see Appendix B), and regardless 
of whether the distribution is gaussian or not. All quantities in the system are related linearly by 
functions of wavenumber fc, because we consider the system at the linear regime. In practice, the 
Fourier transform of the perturbations in each quantity are a time and space dependent amplitude, 
multiplied by a normal Gaussian random number (with or without correlations, depending on the 
type of initial conditions). 

In summary, in the case of adiabatic initial conditions, all quantities in the vector fj:(T) share 
the same random seed, 


%(^) = 




i X. 


■>k ’ 


[adiabatic] 


(3.26) 






regardless of whether these quantities describe a lagrangian displacement field in a mesh which appears 
curved in Eulerian coordinates, or whether these quantities are on a regular Eulerian grid. Already 
in a ACDM universe, baryons and Cold Dark Matter have slightly different transfer functions, which 
should be taken into account when generating their initial conditions based on the same random seed, 
as first applied in [46] . Owing to the linearity of the system, any type of isocurvature perturbations 
can then be expressed as. 


%(^) = 


^ Di (^r, itj ^ 
^ Dm (t, 




^ Di (|r, ^ 

^ Dm j 




[adiabatic + isocurvature] 


(3.27) 


and so on. These statements hold for any type of initial seed, whether it is gaussian or not. Notably, 
in [47] it is pointed out that exactly a scalar modification of gravity may acquire its own spectrum of 
perturbations, giving rise to isocurvature perturbations. 


4 Non-relativistic initial conditions 

4.1 Non-relativistic limit 

The non-relativistic limit for perturbations in an expanding universe is obtained when velocities are 
small compared to the speed of light. Furthermore, the Newtonian limit is obtained when pressure is 
small, in which case the linearised perturbations of the fluid are described by [48], 

h.p = — 9 + q, [Newtonian, ct = 0, w = 0, ^ 1] (4.1) 

0 — q =k^(j)dk^Ap — H{9 — q), (4.2) 

where we continue employing the usual co-expanding® coordinate system with conformal time. How¬ 
ever, these equations neglect the effect of pressure on the energy density present in relativity (associ¬ 
ated to the work needed to compress a fluid with pressure). The equations can be modified to obtain 
the non-relativistic limit that has the correct continuation to the General Relativistic equations [49], 
which in practice amounts to removing all references to geometry and installing a single gravitational 
potential. 


Ap =(1 -I- w){q — 9), [Non-relativistic, cr = 0, re = ^ 0, zi; = 0] (4.3) 

e-q=k^^+^^^k^Ap-H{9-q). (4.4) 

(w -I- 1) 

The Newtonian limit is necessary for setting initial conditions for Newtonian simulations. Power 
spectra for the distribution of densities, velocities and potentials are however readily obtainable from 
solvers of relativistic Boltzmann equations, such as Class [26] and Game [25]. 

One could use the gauge freedom to hx the gauge such that Eqs. (4.3, 4.4) hold at all scales [50, 
51]. We briefly elaborate on that in Section 4.2. In this subsection we focus however on the weak field 
limit, in which initial conditions for Newtonian N-body simulations are inevitably set. An approximate 
translation from Newtonian large-scale results to relativistic interpretations can be found in [52]. 

The proof that Newtonian gravity is the weak field limit of General Relativity is textbook mate¬ 
rial, where it is usually shown that the equations of motion in the longitudinal gauge {Ht = B = 0) 
reproduce the non-relativistic equations of motion (4.3), although the same can be proven in the 
comoving gauges {9 = B) [48]. The density perturbations of both these gauges, computed in terms 

®Nomially called comoving coordinates, we exceptionally say co-expanding in order to avoid confusion with a coor¬ 
dinate gauge for relativistic perturbations, such as the comoving gauge. 
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of quantities of an arbitrary gauge, are [42], 

A^™m°v) + 3(1 + w) ^ (0 - B), (4.5) 

/c 

- . (4.6) 

The pre-factor T-L/k is by definition a measure of the smallness of velocities in the problem^*^, such 
that T-L/k ^ Q gives the lowest order terms in the Newtonian limit. In other words, the density 
perturbation computed in any gauge reduces to the same density perturbation in the non-relativistic 
limit, and is a solution to the set Eqs. (4.3, 4.4), 

Ap = + 0{l) = A^'°"s) + 0{l). (4.7) 

This means that for the non-relativistic limit, one can use the density perturbations computed in any 
gauge (also gauges in which 0 = 0 ), and set the velocities by, 

^(newt) _ ~^P 

1 -b rc ’ 

More generally, this equation holds under any extension of general relativity (modified theory of 
gravity), provided that the fluid in question satisfies Ap = —(1 -|- w)6 with w = 0. 

To summarise: 

• for a newtonian simulation, positions are obtained from the density-displacement duality applied 
to Ap, and velocities from Ap, 

• while for a relativistic simulation, the density-displacement duality is applied to bare quantities, 

A^“® and = -9. 

4.2 Gauges and the Newtonian limit 

From comparing Eqs. (3.14, 3.15) to Eqs. (4.3, 4.4), it is obvious that a gauge in which = 0 
reproduces the same linear equations in General Relativity as in Newtonian gravity. It is tempting to 
believe that Newtonian simulations hence can be used on relativistic distances d such that dTi = 0(1). 
One strong hint that this is the case, comes from the fact that linear perturbation theory agrees well 
with observations of large scale structure in the universe, suggesting that small-scale nonlinearities 
do not affect the large scale dynamics. The growth of structure is hierarchical: small scales enter 
the nonlinear regime first, while large scales continue to evolve linearly. Moreover, all modes follow 
a sequence of (1) relativistic dynamics (super Hubble), (2) newton linear dynamics, (3) nonlinear 
dynamics [33]. While this suggests that long wavelength modes in Newtonian N-body simulations are 
solved for properly, this does not at all mean that the full nonlinear dynamics of the short wavelength 
modes are the same for relativistic and Newtonian systems. What is needed is a proof that the fully 
nonlinear equations of Newtonian gravity and General Relativity agree, even when taking into account 
all scales (since obviously, when any scale goes nonlinear, formally the entire system is nonlinear). It 
is often claimed that even the nonlinear density solutions only source the relativistic potentials at the 
linear level, but as pointed out in [52], it is inconsistent to linearise the left-hand side (the gravity 
part) of the Einstein equations without linearising the right-hand side (the matter part). When the 
matter part contains nonlinear contributions, the average of the matter density will no longer agree 

the Newtonian picture, cosmic expansion is a recession velocity Ureceasion = Hd, for arbitrary distance d. 
Identifying k with a wavelength A = jk, we have Ureceasion = 1 for k = 'H/(27r), hence for simulations that include 
scales close to the Hubble radius. Thus, the Newtonian limit is valid for simulations of boxes smaller than 27r/(3(l -b 
■w)H) [33]. 

'^^During the final stages of preparation of this manuscript, the same was pointed out in Ref. [53], where it is shown 
that for Cold Dark Matter, the gravitational potential follows the Newtonian Poisson equation, such that the linear 
dynamics look completely Newtonian at relativistic scales. 


3(l-bw)H 


< 1 


(4.8) 


Ar'> =A, + 3(1 + »)M^ 
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Figure 2. A 3x5 of patch of a 128^ slice of a realization 128® vertices (particles) of a Dark Matter density 
field at a comoving size of 200 Mpc (1.56 Mpc inter-vertex distance), at newtonian displacements varying from 
redshift 2 = 100 down to 2 = 10, in a cosmology endowed with designer f{R) gravity in wCDM background 
{w = —0.95) with Bo = 0.9, compared to a plain ACDM cosmology. The f{R) trajectories can be distinguished 
from the ACDM counter parts, by their curved shape, as emphasised in the inset. The straight dashed line 
connecting start and end points serves as a guide to the eye, to recognise the curved shape of the trajectory. 
This displacement field is linear, yet the particles do not follow straight lines. 


with the definition of the background matter density, around which perturbation theory was setup. 
Properly taking this mismatch into account, amounts to the Buchert formalism (see Ref. [54] for a 
review), which schematically can be summarised as, 

Gy,u{gt,u{T)) =r^^(p(r), P(r),...), (4.9) 

xj) =5 T^i,{5p{t, x),SP{t, x),. . [ordinary perturbation theory] (4.10) 
^))) ...), [averaging vs. background] (4.11) 

where both the left and right-hand side of the last line are pure functions of r, but different pure 
functions of r. 

In General Relativity, intuition tells that wavelengths larger than the scales considered, can 
be included as a locally constant spatial curvature term (the universe locally is open or closed). 
However, this only addresses the density contributions, while the Einstein equations contain higher 
order gradients in the potential as well. Under general modifications of gravity, Birkhoff’s theorem 
no longer holds, and the above intuition is flawed. 

In summary, it is not obvious that two perturbation theories that are identical in their linear 
limits, are identical in their full nonlinear regimes. 

In [55] is was found numerically that the mismatch between average and background quantities 
may be tiny, but a rigorous mathematical proof of the agreement between fully nonlinear Newtonian 
gravity and General Relativity in cosmic structure formation does currently not exist. 

5 Scale dependent growth in the quasi-static approximation 

In this section, we show that for Dark Energy and Modified Gravity models (in short DE/MG models) 
the velocities of Dark Matter particles are scale dependent in linear perturbation theory, already in 
the quasi-static approximation. Scale dependence of the density transfer function of Dark Matter in 
Eourier space generally translates into a varying direction in real space. That is, the dust particles do 
not move on straight lines. Here we choose to apply the quasi-static approximation in order to make 
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an intuitive transition from General Relativity to the more general Effective Field Theory of Modified 
Gravity. In the next section we abandon the quasi-static approximation^^. 

We fix the metric to the longitudinal gauge, B = Ht = 0, ^ = 'h and Hl = —d>, 

ds^ = a^{T)\—{l+ 2'^)dT^ + {1 — 2^)dx^\ . (5.1) 

Via the quasi-static approximation, we can parametrise the Poisson and anisotropic stress equations 
into the following generic form [56]: 

= AnGa'^ Q{t, x)piAi{T, x) , (5.2) 

-V^(5- - m) = UnGa^p^il + u>i)cr,Q(r, x) , (5.3) 


where Q{t,x) and R{t,x) describe the variation of Newton’s constant in the environment and the 
anisotropic stress tensor induced by the modification of gravity, respectively. The quasi-static ap¬ 
proximation was proposed via {p, 7 ) functions in [57] , but mixing the gravitational shear and the 
anisotropic stress from relativistic species^^. Several phenomenological parameterisations of these 
functions exist in the literature. We refer to the Refs. [58-60] for details. 

Gonsider a system with some parametrisation of modified gravity and a dust fluid (Cold Dark 
Matter). The continuity and momentum equations of the dust fluid in the quasi-static regime (ne¬ 
glecting the time derivatives of the gravitational potentials), can be condensed into. 


d^Ae(T,f) 


+ H 


dAc{T,x) 

dr 




(5.4) 


which can be compared to Eqs. (3.14, 3.15) by setting w = cl = q = Ap = B = a = 0, A = Af and 
$ = 0, where the last equalities follows from the quasi-static approximation. Since, for simplicity, we 
only consider the dust and gravity sectors, the anisotropic stress term from relativistic components, 
such as massive neutrinos, etc. vanishes. 

Combining equation (5.2) and (5.3), and keeping only terms at linear order, we get 


~ —ATTGa^Q{T,x)R{T,x)pcAi.{T,x) . 


(5.5) 


Inserting equation (5.5) into (5.4), we arrive at the equation governing the growth of the dust fluid 
perturbations [61], 


d^Ac(T,x) 


+ n 


dAc{T,x) 

dr 


ATTGa^Q{T,x)R{T,x)pcAc{T,x) = 0 . 


(5.6) 


In General Relativity, Q(t,x) and R{t,x) are unity, hence the perturbations of the dust fluid grow 
in the same way on all the scales, i.e. the growth rate is only a function of time. However, this is no 
longer true in the case of modified gravity. 

The next step is to calculate the growth rate V of CDM density perturbations. In GR, in the 
late-time and small scale, the growth rate only depends on time; while in the modified gravity case, 
'D{t,x) depends both on time and space, ie.. 


Ac(t, x) = D{t, x)Ac{Ti, x) , ( 5 . 7 ) 

where Ac{Ti,x) denotes for the CDM density perturbation at initial time r^. In order to set initial 
conditions for a Newtonian N-body simulation, one now uses Ac for the particle positions. Moreover, 
by ‘virtue’ of the quasi-static approximation, the time derivatives of the potentials are ignored, such 
that bare quantities are equal to the full physical quantities. In other words, both Newtonian and 
relativistic (longitudinal gauge) particle positions are obtained from, 

X = y -V{T,y)^^yAc{Ti,y) - A^{Ti,y)^^yV{T,y) . ( 5 . 8 ) 

our numerical calculation, we do not assume the quasi-static approximation, and evolve the full dynamical 
equation of motion of the extra scalar field. 

^^In the only CDM-|-modified gravity case, the functions (/i-, 7 ) are related to {Q^R) by: Q = /i 7 , R = 7 ”^ . 
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Equation (5.8) is the main result of this section. If we assume that T) depends on time only, we 
recover the Zel’dovich Approximation, which tells us that in General Relativity, at the linear scale, 
the dust particles propagate on a straight line, because {y = Va;Ac(a?); however, in the modified 
gravity scenario, the dust particles are also deflected by the gravitational potential due to the second 
term of equation (5.8). 

Finally, in Figure 2 we show a discrete realisation of Dark Matter particle trajectories for a par¬ 
ticular choice of EFT parameters (see the following section), without the quasi-static approximation. 
The particle trajectories are compared their ACDM counter parts, and clearly show a curved shape. 
This figure shows how particle trajectories at the linear level are deflected by modifications of gravity. 

6 Effective Field Theory parametrisation, and modified gravity as a par¬ 
ticle displacement field 

6.1 Transfer functions 

After having discussed the physical picture, now let us go through a more advanced parametrisation 
method of DE/MG models. In this section, we will focus on the effective field theory (EFT) treatment 
of the cosmic acceleration. This approach is introduced into the study of DE/MG models by [62, 63] 
to unify the most generic viable single scalar field models of DE/MG. Let us briefly summarise it 
here. 

Gompared with the covariant approach, the construction of the effective action of the scalar field 
starts with a particular choice of time coordinate, in which the scalar field perturbations vanish. In 
other words, we first break the four-dimensional covariant diffeomorphisms by choosing a particular 
clock. Then, we can build all the operators that are consistent with the unbroken symmetries of the 
theory, i.e. time-dependent spatial diffeomorphisms. Thanks to this symmetry guidance, we can figure 
out the spatial structure of these operators. Hence, we are left with only time-dependent prefactors 
of these operators, namely EFT functions, which need to be parametrized. Another advantage of this 
procedure is that at the linear order there exist only a few relevant operators which could cover most 
of the generic viable single scalar field DE/MG models. For these reasons, the EFT approach makes 
analysing DE/MG models with the on-going and up-coming cosmological surveys feasible. 

The EFT approach relies on the assumption of the validity of the weak equivalence principle 
which ensures the existence of a metric universally coupled to matter fields and therefore of a well- 
defined Jordan frame. The EFT action in conformal time reads 


S = 


ct^xy/^ I [I -I- 0(t -I- tt)] i? -I- A(r -I- tt) — c(t -I- 7r)a^ 

• 2 

+2TrSg°° -I- 2g°^di'K — ^ + g^’^dindjir — { 2 %^ -|- 'ii^ 


(5g™-22- -p2H7r ( 5g 


„oo 


^21 


Sm Xi] 



( 6 . 1 ) 


where TOq is the Planck mass, and is the action for all matter fields, Xi- For simplicity, here we 
only list the three operators {H,A,c} which is describing the background dynamics. For the complete 
description of the quadratic EFT action, we refer to [27, 62, 63]. In the following calculation, we 
follow the convention of [27]. 

The modified Einstein equation can be written as 

mg(l + L!)G^,[g^,] =r^]7^K,0™,--.]+T^:)[^,7T,...] , (6.2) 

We can divide both sides of the above equation by a factor (1 -|- H) and define an effective energy 
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Figure 3. The fractional differences of the CDM density transfer function between GR and the viable 
f{R) gravity for several redshift snapshots. (Left panel) Designer f{R) gravity in ACDM background with 
Bq = 0.001 and (Right panel) Designer f{R) gravity in wCDM background (w = —0.95) with Bq = 0.01. The 
grey band denotes for the 5% regime. 




Figure 4. The Q-fluid energy density fluctuation transfer function in two viable f{R) gravity. (Left panel) 
Designer f{R) gravity in ACDM background with Bo = 0.001 and (Right panel) Designer f{R) gravity in 
wCDM background (w = —0.95) with Bq = 0.01. 


momentum tensor (T^?^), which is conserved by construction = 0)^^ 

, (6.3) 

TlS'> K, 0., • • • ] = ■ (6.4) 

The reason why we introduce this effective energy momentum tensor (T^^^) instead of directly solve 
the modified Einstein equation is that we want to apply the Lagrangian perturbation treatment for this 
effective DE/MG fluid, hereafter Q-ffuid. The motivation of our algorithms are mainly the following 
two. 

First of all, in an N-body simulation, a high resolution of the extra scalar field in Eulerian 
representation is quite expensive, especially for the collapsing DE/MG models. This is because in 

^^There exists another equivalent definition by moving the non-minimally coupling term from the left-hand side to 
the right, i.e. [pir, Stt, ■ ■ ■ ] = [tt, tt, • ■ ■ ] — [p^^]. The differences between this definition and the 

one in equation (6.4) is that in the former, the extra argument of is a metric field, while in the latter it is a matter 
field. 
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Figure 5. The Q-fluid bare number density fluctuation transfer function of the designer f{R) gravity in 
wCDM background {w = —0.95) with Bo = 0.01. 


the DE/MG simulation, whether we model the extra scalar degree of freedom as the particle or fluid, 
depends on the comparison of the mean free path of the particles with the scales we are interested 
in. Take cold dark matter and massive neutrino as examples, on the scales which are much larger 
than their mean free path, we can adopt the fluid approximation (ideal fluid for CDM; imperfect fluid 
for massive neutrino). However, once we concern on the scales which are comparable or even smaller 
than their mean free path, we have to use the particle description to study its non-linear dynamics. 
A similar case happens to the collapsing DE/MG model, which has a small sound speed (large mean 
free path). One solution is to discretise the fluid element into virtual particles with a charge, e.g. 
mass, and let the grid follow these virtual particles (Q-particles), i.e. the Lagrangian perturbation 
approach. This gives the simulation a high resolution at the high density regions. 

Second of all, the algorithms to add extra fluid components into N-body and hydro simulation 
are being extensively developed. We can utilise these techniques to develop the DE/MG simulations 
via this Q-fluid description, such as loading pressure, etc. In this work, we focus on the linear 
phenomena, such as linear structure growth and the IG for N-body simulation. In principle, this 
Lagrangian treatment can be extended to the non-linear phenomena in the simulations. Basically the 
recipe is to generate the initial conditions for the Q-particles from the linear Boltzmann code; and 
then, evolve them with the geodesics. 

Let us go back to the linear perturbation description. Within this approach the above equations 
(6.3) and (6.4) can be split into two sets of equations, namely background and (linear) perturbation. 
On the background, due to historical reason, pg and Pq are not the conserved background quantities 
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in the non-minimally coupled case, 


Pq = 2c-A -0^ 


The conserved one are defined as 


PQ = A+^{n + nn). 


Pq 


-Pde = — - 


n 


i + n 


P-m + 


P, 


Q 


1 + fl 


(6.5) 

( 6 . 6 ) 


(6.7) 

( 6 . 8 ) 


At the linear perturbation level, equation (6.4) and (6.3) reads 

mlSGf,„[6gfj_„] = 6»„, • • • ] + (5p„, • • • ] , (6.9) 

<5TW) [Sp^, 9^,6p^,---] = [tt, TT, ■ ■ • ]} . (6 .10) 

Armed with these results, we could recognize the fluid variables, such as the energy density, velocity, 
pressure as well as anisotropic stress tensor. Within the fully relativistic treatment, the definition of 
these quantities are gauge related. In the synchronous gauge, the Einstein equation reads 




( 6 . 11 ) 

k fj — {pm + Pm)dL^ ^ + (PDE + Pe)e)Pq ^ , 

( 6 . 12 ) 

-'^{h + 2 nh- 2 k^p^^ 

= 3sp^y’^'> + 36P^y’^^ , 

(6.13) 

A + 67 ) + 2 'H{h + 61 )) — 2 k^p 

= {Pm + Pm) 0 '^^’^^ + {PDE + 7A)e)o’q^”^ , 

(6.14) 


where the fluid variables are defined as 


Sp. 


(syn) 


1 


(1 + fI) 

4 


2 mg 


(pde + Pde)(^q^’^^ = 


1 


1 + fI L 

2 mn 


+ pqtt + 2 c ( 7 T (^ y ") + 

h + j (3(377^ - 

-n(pm + + (pQ + 




+ PQTT^^y^^ + {pq + Pq)(^("^") + 

^ + (ii + 3'Hfi)7T("y") + (nn + + ho + 

r [o \ 3 


a 


{PUE + i^DE)ag""^ = Y^ -^^Pm + + 67 ) + 


(6.15) 


(6.16) 


T^lsyn) 

(6.17) 

(6.18) 


Beware that in the ACDM background (pde + -Pde = 0), the divergence of the velocity and 

the anisotropic stress tensor ctq^”^ of the Q-fiuid are not well defined. The Einstein equation and the 
Q-fiuid energy momentum tensor in the longitudinal gauge are listed in the appendix A. 
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In the rest of this section, we will evolve the full dynamics of the Q-fluid (or tt field) by using 
EFTCamb [27, 28], calculate the relevant fluid quantities for the simulations and produce a realisation 
of discrete initial conditions for a discrete simulation. In the following numerical calculations, for 
simplicity, we take the f{R) gravity as an example to show the non-trlvial dynamics of the Q-fluid. 
In details, we Investigate the designer f{R) model [64, 65] In ACDM and icCDM background {w = 
Pbe/pue), which could exactly reproduce the background history we fixed a priori. After fixing the 
background expansion history, the extra degree of freedom of this higher derivative gravity theory 
still allows for one extra parameter, here we choose It as the present value of the Compton wavelength 
of the scalar field, namely Bq ~ (in Hubble parameter unit) [ 66 ], where fa = Or/. 

We take the Bq parameter value inside the viable regime after Planck-2013 results [28], namely for 
designer f{R) model in ACDM case we take Bq = 0.001; for designer f{R) model In icCDM case we 
take {Bq = 0.01,ru = —0.95). 

In Figure 3 we show the fractional differences between GR and the designer f{R) models, for 
the cold dark matter density perturbations at 11 redshift snapshots. The left panel for the ACDM- 
mimicking case with Bq — 0.001 {w = —1) and the right panel for the wCDM-mimicking case with 
{Bq = 0.01, w = —0.95). The grey band in the figures is the 5%-deviation regime. We can see that 
for the ACDM case at redshift z = 50 (pink curve), where generally N-body simulations start to run, 
the differences are deep inside the 5% regime. However, for the wCDM case on scales smaller than 
(A: > 1 [/i/Mpc]) the differences are outside of the grey band. 

In Figure 4 we show the transfer functions of the energy density of the Q-fluid. Our calculation 
demonstrates that for the ACDM case, dpg is anti-correlated with Spdni on all the scales, i.e. In 
the CDM over-dense region ((Jpcdm > 0) the Q-fluld energy density perturbation keeps under-dense 
{SpQ < 0). This reflects the fact that the CDM perturbation in the f{R) gravity is enhanced in the 
linear regime. To explain this, we take the Poisson equation in the Newtonian limit 

, 2 , 167rG - 6 R 

V” = g <^/Ocdm-^ • (6.19) 

In the CDM over-dense region, iJpcdm > 0, the Newtonian potential ip < 0 and 6 R > 0. In our Q-fluid 
language, equation (6.19) reads 

-k^ip = ^ {Spcdm + Spq) . (6.20) 

Comparing these two equations, we get Spg oc —SR and this quantity stays negative. In order to 
generate the same depth of the gravitational potential well {ip), we need extra CDM fluctuations to 
compensate the negative Spg. This indicates the fact that the growth rate of CDM get enhanced. 
From the modified gravity point of view at the linear scale^^, the effective gravitational constant Geff 
is enhanced by a factor 4/3 compared with those in GR. What we find is nothing new, just another 
explanation of the same phenomena in term of an exotic fluid component. From the right panel of 
Figure 4 we can see that in the deep redshift, the Q-fluid energy density perturbation change the 
sign (on the large scale Is positive; while on the small scale Ap g is negative^®), but in the late 
redshift it does not. This is due to the complicated competition relationship between the Spm and tt 
field In the definition of Spg in equation (6.15). 

In Figure 5 we show the transfer function of the bare number density perturbation (3.20) of the 
Q-fluid (Ajj^^g ). As we have discussed in section (3.4), this quantity defines the positions of vertices at 
rest with the flow of the Q-fluld, in relativistic simulations. From the equation (3.21) we can see that 
Ajj^g® is the integration of Oq over the conformal time. Furthermore, from equation (6.16) we know 
that in the ACDM background 9q is not well defined, neither is Aj^'^g . Given these considerations, 
in Figure 5 we only show Ajj^-g In the rcCDM background case. We can see that, unlike Ap^cdm/fc^, 
the quantity A]]®g/fc^ increases on small scales. This also reflects the fact that the modification of 
gravity in f{R) models becomes more significant on smaller scales, though above the screening scale. 

^^Here the “linear scale”, we mean the scale above the screening scale via the chameleon mechanism [67]. 

^®Our the above explanation is still valid in the niCDM case on the small scale, where the collapsing of Dark Matter 
will happen. 
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Relativistic and Eulerian 
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Q-fluid (Eulerian) 


Newtonian and Lagrangian 



Q-particles (Langrangian) 



Combined 


Figure 6. Slices of realizations of initial conditions at z = 50, for use in relativistic (top row) and newtonian 
(bottom row) simulations, based on linear perturbation solutions from EFTCamb. The setup is identical 
to the /(-R)-cosmology in Figure 2. The space between vertices is 1.56 Mpc. Sizes of the squares represent 
masses of the particles (densities at the vertices), while for the Q-fluid, the color additionally represents 
internal energy (pressure at the vertices), where the colour green indicates a low pressure, which simply traces 
under-densities. We omit the combination ‘Relativistic and Lagrangian’, since CAMB works in synchronous 
gauge, and the velocity of Dark Matter in that gauge is zero, and a Lagrangian representation in that gauge 
is meaningless. The Q-fluid can be treated in the Lagrangian Representation in a relativistic simulation in 
coordinates synchronous and comoving with Dark Matter. A relativistic simulation of Dark Matter only 
makes sense in a gauge that does not produces caustics, such as for example the longitudinal gauge. The top 
row serves as a proof of concept. The displacement of Dark Matter is exaggerated by a factor of 20 in the 
bottom row, for the purpose of contrast in the illustration. 


6.2 Initial conditions for discrete simulations 

In Figure 6 we show an application of all the findings in this paper. We discretise the results of the 
previous sections using a discretisation as described in Appendix B. In the top row, we setup initial 
conditions for a relativistic simulation in synchronous comoving coordinates, as a proof of concept. 
Such a simulation would make little sense, as the motivation for a simulation with Dark Matter is to 
be able to trace the particles beyond shell crossing, which invalidates the synchronous gauge. The 
figures show a regular grid, which is comoving with the Dark Matter. Proper masses (not bare!) are 
indicated by the size of squares. Moreover, for the Q-fluid a low pressure is indicate the colour green, 
while high pressure is indicated in blue. For this fluid, pressure simply traces density. 
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In the bottom row of 6, we show the Newtonian approximation based on the same spectra 
computed in the synchronous comoving gauge, as explained in Section 4. In this case, both fluids can 
be considered in the Lagrangian representation, that is, represented by vertices (particles) which are 
at rest with each species’ flow. From this picture it becomes explicit that the phrase ‘at rest with the 
flow’ does not imply that mass- and internal-energy densities are constant in such a frame, if pressure 
or heat transfer are present. For Dark Matter, pressure is zero and hence these particles do have 
constant mass (equal-sized squares) in the Lagrangian representation^^. 

In the Eulerian view, it is clearly seen that the Q-fluid has density perturbations which are anti¬ 
correlated with the Dark Matter perturbations, since both fluids are over-dense in complementary 
regions. This exemplifies the discussion at Eqs. (6.19, 6.20). However, it is important to realise that 
this regular grid is in the frame comoving with the Dark Matter, which is hence not a regular grid 
in newtonian coordinates; in the newtonian representation (bottom row), one can recognise that the 
low-mass and low-pressure vertices in the Q-fluid are in fact clustered, partially compensating the 
density perturbation that is visible in the top row. 

7 Conclusion 

We have provided the means to generate a distribution of particles with masses, positions and internal 
energies such that they describe exactly the underlying fluid theory, for both perfect and imperfect 
fluids, in arbitrary gauges. That is, simulations in longitudinal, synchronous co-moving, or the so- 
called ‘N-body gauge’ [53], find their initial conditions following our prescription. The newtonian limit 
is always found by taking the dressed densities for particle positions, while fully relativistic particle 
positions are defined by the bare number densities. The description holds for any metric theory of 
gravity, including the newtonian limit of General Relativity. The description can be used in any 
gauge, for any description of a fluid, with velocities defined as comoving with any quantity of the 
fluid. 

We applied our findings to a description of Modified Gravity, in which we describe the extra 
degree of freedom as an imperfect fluid Q, with particles with varying mass and internal energy. 
Trajectories of all species in the linear regime are curves rather than straight lines. Most importantly, 
we show that the perturbations in new dynamical degrees of freedom cannot be ignored when setting 
initial conditions. 

Such initial conditions can be applied in at least three ways: (I) when the modified gravity field 
stays linear and its fluid description remains valid, a Lagrangian representation of the linear field can 
be included in an N-body simulation during all times, providing high resolution where the simulation 
needs it, (2) when the nonlinear dynamics of the field are still described by a fluid, the field can be 
simulated by adapting existing hydrodynamics techniques, and (3) these initial conditions can be used 
for any particle species whose linear solution is effectively described by an imperfect fluid, such as 
neutrinos on scales larger than their free-streaming length. 

Our findings open the way to modelling non-Newtonian gravity in various ways, in arbitrary 
coordinate choices. The advantages of this possibility remain to be explored. 

We release a G+^-code for the generation of initial conditions, EalconIC, at http: //f alconb. org, 
with a minimalistic GUI that runs natively on Linux and OS X. The code links against any version 
of both Boltzmann codes CAMB and CLASS, including EETCamb, such that no separate running of 
those codes is necessary, and generates initial conditions at arbitrary scales, of arbitrary size (fully 
parallelised using MPI and OpenMP), for arbitrary cosmological parameters. 

Acknowledgements 

The authors wish to thank Ignacy Sawicki, Julien Bel, Garmelita Garbone, Julien Lesgourgues, Mark 
Lovell, Ewald Puchwein, Gornelius Rampf, Marco Raveri, Gerasimos Rigopoulos, Christoph Schmid 
and Alessandra Silvestri for fruitful discussions. WV is supported by a Veni research grant from the 

there is a direct coupling between Dark Matter and the extra scalar field, the Dark Matter might not be conserved 
due to the non-trivial coupling[68— 71]. 


- 19 - 



Netherlands Organization for Scientific Research (NWO). BH is supported by the Dutch Foundation 
for Fundamental Research on Matter (FOM). 

A Gauge transformation and Q-fluid variables in the longitudinal gauge 

In the longitudinal gauge, with the convention of Ma and Bertschinger [40] , 

ds^ = ~ (1 + ‘^^)dT'^ + (1 ~ 2(j))5ijdx^dx^ 


ds^ = a? — (1 + 2A)dT‘^ + (1 + 2H]^)5ijdx'^dx^ 


so we have A = ip , = —</>■ In the following, we will move to the {^p, p) convention 

1 




h + Qrj + 'H{h + 677 ) 
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^(long) ^ 

In the longitudinal gauge, the Einstein equation reads 
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And the Q-Auid variables are defined as 
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B Discrete sampling of phase space 


For the sake of completeness, here we briefly summarise the currently known methods for generating 
a discrete sampling of the otherwise continuous phase space of a fluid. For in-depth discussions, we 
refer the reader to [45, 72], 

In this section we work top down, in that we start from a theoretical continuous field on an 
infinite space, and construct a finite-resolution finite-sized discrete representation that resembles the 
theoretical field as closely as possible. 

Discreteness enters at two ends of the length scale: (1) moving from an infinite box size to a finite 
box size renders the eigen space of the laplacian discrete (Fourier transforms become sums rather than 
integrals), which sets the lowest nonzero wave number under consideration; and (2) only finitely many 
numbers can be treated on a discrete computer, which on one hand limits the highest sampled wave 
number and on the other hand limits the number of simulated vertices (particles). 

A statistically isotropic gaussian random field s(x) = f in an infinite-sized space can 

be defined by its Fourier-space correlator. 


(%%,) = (27r)353(fc-fcV^(fc) 


(B.l) 


such that 



(B.2) 


(B.3) 


is independent of position x. Random variables are denoted by a hat, s. 


B.l Finite length 

The continuous field s(x) in the hnite-sized box with length L and periodic boundary conditions (a 
3-torus) is related to the continuous field in infinite-space. 



OO 


= lim 

L—^oo 



(B.4) 


where n is a vector of three integers. Inserting Eq. (B.4) into Eq.(B.2), one finds that. 



(B.5) 


where is the Kronecker delta. 


B.2 Finite resolution 

Next, keeping the finite box size L (not taking the infinite limit), the relation between the continuous 
and discrete functions in the box is set by the highest wave number ttM/L, 


M/2 


s(x) = lim 

M-foo 



(B.6) 


n=-M/2 


It is custom to simply take a finite M, which is equivalent to multiplying the Fourier space Sa^rs 
with a tophat filter. Different filters may be chosen. The real space representation of the tophat 
filter indeed closely resembles a series of unit valued peaks, equally spaced at a separation L/M, such 









that the finite M leads to a real space s{x) which is convolved with a series of delta functions, i.e. 
a discrete sampling. Note however that relation between the tophat filter and a discrete real space 
sampling is not exact, though it is a good approximation. Formally, a discrete periodic sampling of a 
function in real space corresponds to a discrete periodic sampling in Fourier space exactly, only when 
both samples are infinite (M —)• oo, see [73]). 

When one wants to generate a coarse grained realisation that represents an averaging of the 
continuous field, the filter applied inside the sum of Eq. (B.6) should be the Fourier transform of the 
filter used in the real space averagin, such as a gaussian filter, cloud-in-cell, triangular-shaped-cloud, 
and so on [74]. 

In summary, the continuous gaussian random field s{x) on an infinite space, can be discretely 
sampled on a grid with length L and resolution M by generating gaussian random numbers that obey. 


{s{nL / M)s{n' L / M)) 


M/2 

rh— — Ml2 
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with n a vector of three integers in the range [0,M) and m a vector with integers in the range 

[-M/2,M/2). 


B.3 Generation in Fourier space 

From Eq. (B.5), we find that both the real and imaginary parts of S 2 jva are independent gaussian 
random variables (see e.g. [14]) with variance Ug = 

2 \ 
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Imposing Hermitian symmetry Sa^ri = s_ 22 ^ guarantees that s(x) be real valued. Einally, one simply 

~L~ i- 

has to draw normal gaussian random complex values ^ for each discrete wave number, and multiply 
by its respective amplitude to realise one random field, 

M/2 


Kxi) = 

n=-M/2 




27rn 




(?) 


(B.9) 


where Xi represents a position on the discrete grid, with three discrete components in the range [0, M). 


B.4 Generation in real space 

Eollowing [45], the separation of amplitude and normal random variables in Eq. (B.9) implies that 
s{xi) can be written as a convolution of the inverse Fourier transforms of the individual ingredients, 
the amplitude and the normal random field. A normal random field transforms into a normal random 
field. Hence, one could just as well generate the normal random field in real space, and perform the 
convolution with the real space transform. 


r(f 0 


M/2 

E ' 

ri=-M/2 


Xi 

\L-^Ps i 

27rn 

)] 


2 SI 

L 

)\ 


(B.IO) 


Generating just one grid, this method is equivalent to the Fourier space method. However, 
now having access to the random numbers in real space, this opens the way to generating higher 
resolution initial conditions in only a part of the initial conditions, by creating a subgrid in a subregion, 
which is again convolved with the real space transfer function T{xi). These are so-called zoom initial 
conditions. The description here is by no means sufficient to generate accurate zoom initial conditions, 
but is intended as a rough sketch. Refs. [45, 72] discuss the matter in great detail, solving issues such 
as mass conservation and keeping the zoom region ‘invisible’ to the coarsely sampled outer regions 
when it comes to the gravitational Poisson equation. 
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B.5 Pre-initial conditions, or particles 

The choice of largest wave number k = ttM/L is unrelated to the choice of the number of particles, 
although it is tempting (and simplest) to set them equal, such that the vertices on the grid for the 
Fourier transform are directly interpreted as particles. This simplest choice is what we used in the 
illustrations in the body of this paper. An obvious deviation is to choose less particles, still on the 
vertices though, simply skipping a fraction of the generated vertices. 

The most popular choice is to put particles on glass pre-initial conditions. Particles are dis¬ 
tributed randomly and evolved under a repulsive force with a friction to aid convergence, until they 
reach an equilibrium. Such a configuration has no preferred directions, as opposed to a grid. The 
downside is that the particles are not at the vertices of the Fourier grid, such that the values of the 
field at the particle positions must be obtained by interpolations, and a mathematically rigorous con¬ 
struction of the held values at particle positions is lost. The interpolation between the grid vertices, 
boils down to turning the discrete held back into a continuous one, by convolving the discrete held 
with a mass window function, the shape of which needs to be taken into account as it affects the 
power spectrum just like a hlter. Other tilings, apart from glass and grid, exist as well [36]. 

As pointed out in for example [12], both grids and glasses are gravitationally unstable against 
small perturbations. This implies that part of the hnest structure found in a simulation is due to the 
choice of pre-initial conditions, and not cosmological of origin. 

B.6 Summary 

The parameters that a simulator has to choose in the initial conditions, include at least: 

• the box size L, 

• the highest wavenumber ■kM/L 

• additional hlters on the initial spectrum (apart from the tophat at ttM/L), 

• the number of vertices / particles, 

• the position of the vertices (on the Fourier grid, glass, ...). 

Each of these has its own consequences for the resolution of the simulation, and the meaning of its 
outcome. What these consequences are, is a whole field of research [11-18], and is beyond the scope 
of this paper. 
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